
options(scipen=20)

library(MASS)
library(WhatIf)
library(RColorBrewer)
library(simcf)
library(verification)
library(foreign)
library(nnet)
library(xtable)
library(apsrtable)
library(ggplot2)
library(dplyr)
library(stargazer)
library(haven)
col <- brewer.pal(3,"Dark2")

beat_data <- read_dta("beat_counts.dta") %>% na.omit()

# basic analysis (Models found in Table 3 of the manuscript) #
m1 <- lm(to14 ~ lcon+lcbo+pct_young+pct_old+pct_black+pct_latino+pct_college+pct_pov+pct_unemp+pct_owner, data = beat_data)
m2 <- lm(to15 ~ lcon+lcbo+pct_young+pct_old+pct_black+pct_latino+pct_college+pct_pov+pct_unemp+pct_owner, data = beat_data)

m3 <- lm(mtg_p1k ~ lcon+lcbo+pct_young+pct_old+pct_black+pct_latino+pct_college+pct_pov+pct_unemp+pct_owner, data = beat_data)
m4 <- lm(call311_p1k ~ lcon+lcbo+pct_young+pct_old+pct_black+pct_latino+pct_college+pct_pov+pct_unemp+pct_owner, data = beat_data)

stargazer(m1,m2,m3,m4, no.space=TRUE,type="text",
          single.row = TRUE,
          covariate.labels=c("Log(Convictions per 1000 Pop)",
                             "log(CSO per 1000 Pop)","% between 18-34","% 65+",
                             "% Black","% Latino", "% College Grad","% Living in Poverty",
                             "% Unemployed","% Owner Occupied",
                             "Constant"),
          dep.var.labels.include = FALSE,
          model.numbers=FALSE,
          star.char=c("*","**","***"),
          star.cutoffs = c(0.05,0.01,0.001),
          column.labels = c("2014 Turnout","2015 Turnout",
                            "Meeting Attendence per 1000 Pop","311 Calls Per 1000 Pop"),
          dep.var.labels   = "")

# moderating effect of CSO density (found in Table 4 of manuscript)  #

m1 <- lm(to14 ~ lcon*lcbo+pct_young+pct_old+pct_black+pct_latino+pct_college+pct_pov+pct_unemp+pct_owner, data = beat_data)
m2 <- lm(to15 ~ lcon*lcbo+pct_young+pct_old+pct_black+pct_latino+pct_college+pct_pov+pct_unemp+pct_owner, data = beat_data)

m3 <- lm(mtg_p1k ~ lcon*lcbo+pct_young+pct_old+pct_black+pct_latino+pct_college+pct_pov+pct_unemp+pct_owner, data = beat_data)
m4 <- lm(call311_p1k ~ lcon*lcbo+pct_young+pct_old+pct_black+pct_latino+pct_college+pct_pov+pct_unemp+pct_owner, data = beat_data)

stargazer(m1,m2,m3,m4, no.space=TRUE,type="text",
          single.row = TRUE,
          covariate.labels=c("Log(Convictions per 1000 Pop)",
                             "log(CSO per 1000 Pop)","% between 18-34","% 65+",
                             "% Black","% Latino", "% College Grad","% Living in Poverty",
                             "% Unemployed","% Owner Occupied",
                             "Convict*CSOs",
                             "Constant"),
          dep.var.labels.include = FALSE,
          model.numbers=FALSE,
          star.char=c("*","**","***"),
          star.cutoffs = c(0.05,0.01,0.001),
          column.labels = c("2014 Turnout","2015 Turnout",
                            "Meeting Attendence per 1000 Pop","311 Calls Per 1000 Pop"),
          dep.var.labels   = "")

## Visualize (Figure 3 of manuscript) ##

coefs = as.data.frame(summary(m1)$coefficients[c(2,3,12),1:2])
names(coefs)[2]="se"
coefs$vars = rownames(coefs)

variable<-c("01","02","03")
d<-cbind(variable,coefs)
p1<-ggplot(data = d, aes(x = variable, y = Estimate, 
                         ymin=(Estimate - 1.96*se), ymax = (Estimate + 1.96*se))) +
  geom_pointrange(size=.5)+
  theme_bw(base_size=16,base_family="Times")+
  scale_x_discrete(breaks=c("01","02","03"),labels=c("Log(Convict)","Log(CSO)","Convict*CSO"))+
  geom_hline(yintercept = 0, lty = 2)+
  theme(plot.margin=unit(c(1,2,1,1),"cm"))+
  labs(y="Marginal Effect\n",x="\n Model: 2014 Turnout",title=c(" ", cex=1))

ggsave(p1, filename = paste("coeff_14TO", ".pdf", sep = ""),width = 7, height = 5)

# 15 turnout #
coefs = as.data.frame(summary(m2)$coefficients[c(2,3,12),1:2])
names(coefs)[2]="se"
coefs$vars = rownames(coefs)

variable<-c("01","02","03")
d<-cbind(variable,coefs)
p1<-ggplot(data = d, aes(x = variable, y = Estimate, 
                         ymin=(Estimate - 1.96*se), ymax = (Estimate + 1.96*se))) +
  geom_pointrange(size=.5)+
  theme_bw(base_size=16,base_family="Times")+
  scale_x_discrete(breaks=c("01","02","03"),labels=c("Log(Convict)","Log(CSO)","Convict*CSO"))+
  geom_hline(yintercept = 0, lty = 2)+
  theme(plot.margin=unit(c(1,2,1,1),"cm"))+
  labs(y="Marginal Effect\n",x="\n Model: 2015 Turnout",title=c(" ", cex=1))

ggsave(p1, filename = paste("coeff_15TO", ".pdf", sep = ""),width = 7, height = 5)

# meetings #
coefs = as.data.frame(summary(m3)$coefficients[c(2,3,12),1:2])
names(coefs)[2]="se"
coefs$vars = rownames(coefs)

variable<-c("01","02","03")
d<-cbind(variable,coefs)
p1<-ggplot(data = d, aes(x = variable, y = Estimate, 
                         ymin=(Estimate - 1.96*se), ymax = (Estimate + 1.96*se))) +
  geom_pointrange(size=.5)+
  theme_bw(base_size=16,base_family="Times")+
  scale_x_discrete(breaks=c("01","02","03"),labels=c("Log(Convict)","Log(CSO)","Convict*CSO"))+
  geom_hline(yintercept = 0, lty = 2)+
  theme(plot.margin=unit(c(1,2,1,1),"cm"))+
  labs(y="Marginal Effect\n",x="\n Model: Attendance at Police Beat Meetings",title=c(" ", cex=1))

ggsave(p1, filename = paste("coeff_meetings", ".pdf", sep = ""),width = 7, height = 5)


# 311 Calls #
coefs = as.data.frame(summary(m4)$coefficients[c(2,3,12),1:2])
names(coefs)[2]="se"
coefs$vars = rownames(coefs)

variable<-c("01","02","03")
d<-cbind(variable,coefs)
p1<-ggplot(data = d, aes(x = variable, y = Estimate, 
                         ymin=(Estimate - 1.96*se), ymax = (Estimate + 1.96*se))) +
  geom_pointrange(size=.5)+
  theme_bw(base_size=16,base_family="Times")+
  scale_x_discrete(breaks=c("01","02","03"),labels=c("Log(Convict)","Log(CSO)","Convict*CSO"))+
  geom_hline(yintercept = 0, lty = 2)+
  theme(plot.margin=unit(c(1,2,1,1),"cm"))+
  labs(y="Marginal Effect\n",x="\n Model: 311 Calls",title=c(" ", cex=1))

ggsave(p1, filename = paste("coeff_311calls", ".pdf", sep = ""),width = 7, height = 5)




